This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.

Attach the required R packages

#library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
## Warning: package 'ggplot2' was built under R version 3.1.3
library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)
## Warning: package 'car' was built under R version 3.1.2
library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)

VIS analysis

Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)

#if (!file.exists('vis_snapshots.txt')) {
#  download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
#                'vis_snapshots.txt')
#}

Read data and format for analysis

vis.data = read.table(file="vis_snapshots_nocorrect.csv", sep=",", header=TRUE)

Planting date

planting_date = as.POSIXct("2013-11-26")

Add water treatment column coded in barcodes.

vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66

Add plant genotype column coded in barcodes.

vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'

Add genotype x treatment group column

vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')

Add calendar-time data column using the Unix-time data

vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")

Calculate days after planting from planting data

vis.data$dap = as.numeric(vis.data$date - planting_date)

Soil volume water content

Use linear regression to create a simple model for using water volume to predict soil volume water content.

Download data for soil wet/dry weight and volume water content measurements.

if (!file.exists('soil_weigth_vwc.txt')) {
  download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
                'soil_weigth_vwc.txt')
}

Read the soil volume water content data

vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)

Calculate the average soil dry weight

mean(vwc.data$weight_dry)
## [1] 72.92138

Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.

vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)
## 
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol < 
##     260, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5426 -1.6196  0.3976  1.4214  4.0771 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.032863   1.389638  -2.182   0.0391 *  
## water_vol    0.233197   0.008145  28.630   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared:  0.9716, Adjusted R-squared:  0.9704 
## F-statistic: 819.7 on 1 and 24 DF,  p-value: < 2.2e-16

Plot the model

vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
                  geom_point(size=2) +
                  geom_smooth(method='lm', formula=y~x) +
                  scale_x_continuous("Water volume (ml)") +
                  scale_y_continuous("Volume water content (%)") +
                  theme_bw() +
                  theme(axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold"))
print(vwc.plot)

Predict volume water contents for the water treatment groups

treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)
##   water_vol      vwc
## 1     217.0 47.57084
## 2     144.5 30.66407
## 3      72.0 13.75731

Convert VIS camera zoom units

LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.

zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
                                                    zoom.camera=c(1,6)))
summary(zoom.lm)
## 
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1, 
##     6000), zoom.camera = c(1, 6)))
## 
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665         NA      NA       NA
## zoom        0.0008335         NA      NA       NA
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 1 and 0 DF,  p-value: NA

VIS zoom correction

In this section we define models that are used to convert area and length between camera zoom levels to a common scale.

Download data for a reference object imaged at different zoom levels.

#if (!file.exists('zoom_calibration_data.txt')) {
#  download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
#                'zoom_calibration_data.txt')
#}

Read zoom calibrartion data

z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)

Calculate px per cm

z.data$px_cm = z.data$length_px / z.data$length_cm

Convert LemnaTec zoom units to camera zoom units

z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)
vis.data$zoom = vis.data$sv_zoom
vis.data$sv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)
vis.data$zoom = vis.data$tv_zoom
vis.data$tv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)

Zoom correction for area

Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)

Non-linear regression (exponential)

area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
                     z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
               data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)
## 
## Formula: rel_area ~ a * exp(b * zoom.camera)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  0.34409    0.01634   21.06   <2e-16 ***
## b  0.88510    0.01206   73.40   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2738 on 33 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 1.513e-08

Non-linear regression (polynomial)

area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)
## 
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7458 -0.4521  0.1182  0.3087  1.6483 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.85322    0.54151   8.962 3.08e-10 ***
## zoom.camera      -5.11648    0.45312 -11.292 1.07e-12 ***
## I(zoom.camera^2)  1.73120    0.08478  20.419  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.496 on 32 degrees of freedom
## Multiple R-squared:  0.9892, Adjusted R-squared:  0.9885 
## F-statistic:  1465 on 2 and 32 DF,  p-value: < 2.2e-16

AIC

AIC(area.nls, area.pol)
##          df      AIC
## area.nls  3 12.58188
## area.pol  4 55.10572

The exponential model minimizes AIC. Plot exponential model.

nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
                  geom_point(size=2.5) +
                  scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
                  scale_y_continuous(lim=c(0,17),
                                     "Reference object relative pixel area") +
                  stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
                              method="nls", se=FALSE, formula=y ~ a * exp(b * x),
                              start=c(a=area.coef$a, b=area.coef$b)) +
                  theme_bw() +
                  theme(axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold"))
nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)

Zoom correction for length

Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).

Non-linear regression (exponential)

len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
                    z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
              data = z.data[z.data$camera == 'VIS SV',],
              start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)
## 
## Formula: px_cm ~ a * exp(b * zoom.camera)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 16.580319   0.271317   61.11   <2e-16 ***
## b  0.449384   0.004676   96.10   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.035 on 15 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 1.958e-06

Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.

len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
              data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)
## 
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera == 
##     "VIS SV", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75431 -0.44016 -0.09013  0.40661  1.28788 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       31.1697     0.9308   33.49 9.13e-15 ***
## zoom.camera       -8.9629     0.7901  -11.34 1.92e-08 ***
## I(zoom.camera^2)   6.5778     0.1503   43.75 2.24e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5886 on 14 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9995 
## F-statistic: 1.718e+04 on 2 and 14 DF,  p-value: < 2.2e-16

AIC

AIC(len.nls, len.poly)
##          df      AIC
## len.nls   3 53.28426
## len.poly  4 34.92045

The polynomial model minimizes AIC. Plot polynomial model.

poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
                   geom_point(size=2.5) +
                   scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
                   scale_y_continuous(lim=c(0,150),
                                      "Reference object length scale (px/cm)") +
                   stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
                               aes(x=zoom.camera, y=px_cm), method="lm",
                               formula=y ~ x + I(x^2)) +
                   theme_bw() +
                   theme(axis.title.x=element_text(face="bold"),
                         axis.title.y=element_text(face="bold"))
poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)

Correct VIS data for differences in camera zoom

Create zoom-corrected VIS data frame

vis.data.zoom = vis.data[,c('plant_id', 'datetime', 'treatment', 'genotype', 'group', 'date', 'dap', 'solidity', 'outlier')]

Predict relative area zoom correction factors

vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$sv_rel_area = predict(object = area.nls, newdata = vis.data)
vis.data$zoom.camera = vis.data$tv.zoom.camera
vis.data$tv_rel_area = predict(object = area.nls, newdata = vis.data)

Calculate total zoom-corrected side-view and top-view area

vis.data.zoom$sv_area = (vis.data$sv0_area / vis.data$sv_rel_area) + (vis.data$sv90_area / vis.data$sv_rel_area) +
                        (vis.data$sv180_area / vis.data$sv_rel_area) + (vis.data$sv270_area / vis.data$sv_rel_area)

vis.data.zoom$tv_area = vis.data$tv_area / vis.data$tv_rel_area

Calculate zoom-corrected lengths

vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$px_cm = predict(object = len.poly, newdata=vis.data)
vis.data.zoom$extent_x = vis.data$extent_x / vis.data$px_cm
vis.data.zoom$extent_y = vis.data$extent_y / vis.data$px_cm
vis.data.zoom$height_above_bound = vis.data$height_above_bound / vis.data$px_cm

Height modeling

In this section we measure how well PlantCV estimates plant height.

Download data manually measured plant height data set (n = 173).

#if (!file.exists('height_model_data.txt')) {
#  download.file('http://files.figshare.com/1845361/height_model_data.txt',
#                'height_model_data.txt')
#}

Read height data.

manual.ht.data = read.table(file="manual_height_samples.csv",sep=",",header=TRUE)

Convert human-readable date to Unix time

manual.ht.data$datetime = as.numeric(ymd_hms(manual.ht.data$timestamp, tz="America/Chicago"))

Get height data from the VIS data for each manual height sample

ht.data = merge(manual.ht.data, vis.data.zoom, by = c('plant_id', 'datetime'))

Convert LemnaTec zoom units to camera zoom units

ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)

Height linear model

height.model = lm(manual_height_cm~height_above_bound, ht.data)
summary(height.model)
## 
## Call:
## lm(formula = manual_height_cm ~ height_above_bound, data = ht.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.99465 -0.23346  0.01155  0.24932  2.21212 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.372426   0.080998  -4.598  7.7e-06 ***
## height_above_bound  1.024824   0.002947 347.694  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6406 on 193 degrees of freedom
## Multiple R-squared:  0.9984, Adjusted R-squared:  0.9984 
## F-statistic: 1.209e+05 on 1 and 193 DF,  p-value: < 2.2e-16

Plot the height linear model.

hlm.plot = ggplot(ht.data, aes(x=height_above_bound, y=manual_height_cm)) +
                  geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
                  geom_smooth(method="lm", formula=y~x, color='black') +
                  scale_x_continuous(lim=c(0,60), "Estimated height (cm)") +
                  scale_y_continuous(lim=c(0,60), "Manually measured height (cm)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.75),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  labs(color="Zoom")
hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)

Biomass modeling

In this section we model fresh- and dry-weight above ground biomass using image measurements.

Download data manually measured plant biomass data set (n = 41).

#if (!file.exists('biomass_model_data.txt')) {
#  download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
#                'biomass_model_data.txt')
#}

Read biomass data.

manual.st.data = read.table(file='manual_biomass_samples.csv', sep=",", header=TRUE, stringsAsFactors=FALSE)

Get data from the VIS data for each manual biomass sample

st.data = merge(manual.st.data, vis.data.zoom, by = c('plant_id', 'datetime'))

Create out-of-frame indicator variable.

st.data$outind = NA
st.data[st.data$outlier == 'True',]$outind = 1
st.data[st.data$outlier == 'False',]$outind = 0

Genotype indicator variable.

st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1

Fresh-weight biomass

Full model. Includes side-view area, top-view area and height.

fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)

Step-wise model selection with AIC.

fw.step = stepAIC(fw.full, direction="both")
## Start:  AIC=-54.68
## fresh_weight ~ sv_area * tv_area * height_above_bound
## 
##                                      Df Sum of Sq    RSS     AIC
## - sv_area:tv_area:height_above_bound  1  0.031061 7.3442 -56.506
## <none>                                            7.3132 -54.680
## 
## Step:  AIC=-56.51
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area + 
##     sv_area:height_above_bound + tv_area:height_above_bound
## 
##                                      Df Sum of Sq    RSS     AIC
## - sv_area:tv_area                     1  0.007503 7.3517 -58.464
## - tv_area:height_above_bound          1  0.059093 7.4033 -58.177
## - sv_area:height_above_bound          1  0.143762 7.4880 -57.711
## <none>                                            7.3442 -56.506
## + sv_area:tv_area:height_above_bound  1  0.031061 7.3132 -54.680
## 
## Step:  AIC=-58.46
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound + 
##     tv_area:height_above_bound
## 
##                              Df Sum of Sq    RSS     AIC
## - tv_area:height_above_bound  1   0.23599 7.5877 -59.169
## - sv_area:height_above_bound  1   0.31978 7.6715 -58.718
## <none>                                    7.3517 -58.464
## + sv_area:tv_area             1   0.00750 7.3442 -56.506
## 
## Step:  AIC=-59.17
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
## 
##                              Df Sum of Sq    RSS     AIC
## - sv_area:height_above_bound  1   0.30538 7.8931 -59.551
## <none>                                    7.5877 -59.169
## + tv_area:height_above_bound  1   0.23599 7.3517 -58.464
## - tv_area                     1   0.56816 8.1559 -58.208
## + sv_area:tv_area             1   0.18440 7.4033 -58.177
## 
## Step:  AIC=-59.55
## fresh_weight ~ sv_area + tv_area + height_above_bound
## 
##                              Df Sum of Sq    RSS     AIC
## - tv_area                     1    0.2631  8.156 -60.207
## <none>                                     7.893 -59.551
## + sv_area:height_above_bound  1    0.3054  7.588 -59.169
## + tv_area:height_above_bound  1    0.2216  7.671 -58.718
## + sv_area:tv_area             1    0.0478  7.845 -57.800
## - height_above_bound          1    0.8110  8.704 -57.541
## - sv_area                     1   23.8156 31.709  -4.536
## 
## Step:  AIC=-60.21
## fresh_weight ~ sv_area + height_above_bound
## 
##                              Df Sum of Sq     RSS     AIC
## <none>                                      8.156 -60.207
## + tv_area                     1     0.263   7.893 -59.551
## - height_above_bound          1     0.654   8.810 -59.046
## + sv_area:height_above_bound  1     0.000   8.156 -58.208
## - sv_area                     1   157.554 165.710  61.263
summary(fw.step)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + height_above_bound, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95619 -0.19476 -0.00613  0.06206  1.49454 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.455e-02  2.002e-01   0.173    0.864    
## sv_area             3.972e-05  1.466e-06  27.093   <2e-16 ***
## height_above_bound -4.029e-02  2.309e-02  -1.745    0.089 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4633 on 38 degrees of freedom
## Multiple R-squared:  0.9842, Adjusted R-squared:  0.9833 
## F-statistic:  1181 on 2 and 38 DF,  p-value: < 2.2e-16

AIC model

fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
              sv_area*height_above_bound, st.data)
summary(fw.aic)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound + 
##     sv_area * height_above_bound, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.14016 -0.10203 -0.00920  0.07919  1.56414 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.812e-01  2.356e-01   0.769   0.4468    
## sv_area                     4.538e-05  4.288e-06  10.584 1.34e-12 ***
## tv_area                    -1.026e-05  6.248e-06  -1.642   0.1093    
## height_above_bound         -6.232e-02  2.709e-02  -2.301   0.0273 *  
## sv_area:height_above_bound  1.924e-07  1.598e-07   1.204   0.2366    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4591 on 36 degrees of freedom
## Multiple R-squared:  0.9853, Adjusted R-squared:  0.9836 
## F-statistic:   602 on 4 and 36 DF,  p-value: < 2.2e-16

The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.

fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10964 -0.16717  0.00494  0.24429  1.30677 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.704e-01  1.002e-01   -2.70   0.0102 *  
## sv_area      3.755e-05  7.931e-07   47.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared:  0.9829, Adjusted R-squared:  0.9825 
## F-statistic:  2241 on 1 and 39 DF,  p-value: < 2.2e-16

Goodness of fit.

anova(fw.aic, fw.red)
## Analysis of Variance Table
## 
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area * 
##     height_above_bound
## Model 2: fresh_weight ~ sv_area
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     36 7.5877                           
## 2     39 8.8099 -3   -1.2222 1.9329 0.1417

SV area model.

sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10964 -0.16717  0.00494  0.24429  1.30677 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.704e-01  1.002e-01   -2.70   0.0102 *  
## sv_area      3.755e-05  7.931e-07   47.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared:  0.9829, Adjusted R-squared:  0.9825 
## F-statistic:  2241 on 1 and 39 DF,  p-value: < 2.2e-16

Plot SV model

sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
                       geom_smooth(method="lm", color="black", formula = y ~ x) +
                       geom_point(size=2.5) +
                       scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                       scale_y_continuous("Fresh-weight biomass (g)") +
                       theme_bw() +
                       theme(axis.title.x=element_text(face="bold"),
                             axis.title.y=element_text(face="bold"))
sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)

SV area model with out-of-frame.

sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)
## Analysis of Variance Table
## 
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1     39 8.8099                                
## 2     38 7.0902  1    1.7196 9.2164 0.004316 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sv.ind.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78369 -0.19699 -0.01344  0.21652  1.20823 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.202e-01  9.252e-02  -2.381  0.02241 *  
## sv_area      3.831e-05  7.636e-07  50.174  < 2e-16 ***
## outind      -5.242e-01  1.727e-01  -3.036  0.00432 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.432 on 38 degrees of freedom
## Multiple R-squared:  0.9862, Adjusted R-squared:  0.9855 
## F-statistic:  1361 on 2 and 38 DF,  p-value: < 2.2e-16

Plot SV model with out-of-frame.

sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight, 
                                        group=outlier, color=outlier)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Fresh-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Out-of-frame")
print(sv.model.out.plot)

SV area model with genotype

sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24908 -0.29687  0.01704  0.21582  1.17723 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.150e-01  1.266e-01  -3.279  0.00223 ** 
## sv_area      3.775e-05  7.797e-07  48.414  < 2e-16 ***
## group        2.614e-01  1.460e-01   1.790  0.08135 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4624 on 38 degrees of freedom
## Multiple R-squared:  0.9842, Adjusted R-squared:  0.9834 
## F-statistic:  1186 on 2 and 38 DF,  p-value: < 2.2e-16

Plot SV model with genotype.

sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight, 
                                        group=genotype, color=genotype)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Fresh-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Genotype")
print(sv.model.gen.plot)

Dry-weight biomass

SV area model

dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16009 -0.04743  0.01569  0.04033  0.18561 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.915e-02  1.389e-02  -3.539  0.00106 ** 
## sv_area      4.425e-06  1.100e-07  40.233  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06591 on 39 degrees of freedom
## Multiple R-squared:  0.9765, Adjusted R-squared:  0.9759 
## F-statistic:  1619 on 1 and 39 DF,  p-value: < 2.2e-16

Plot dry-weight side-view model

dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
                           geom_smooth(method="lm", color="black", formula = y ~ x) +
                           geom_point(size=2.5) +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold"))
dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)

SV area model with out-of-frame.

dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.112715 -0.030983  0.006637  0.024625  0.127446 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.052e-02  1.196e-02  -3.388 0.001650 ** 
## sv_area      4.556e-06  9.869e-08  46.166  < 2e-16 ***
## outind      -9.023e-02  2.232e-02  -4.043 0.000248 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05583 on 38 degrees of freedom
## Multiple R-squared:  0.9836, Adjusted R-squared:  0.9827 
## F-statistic:  1136 on 2 and 38 DF,  p-value: < 2.2e-16
dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight, 
                                        group=outlier, color=outlier)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Out-of-frame")
print(dry.sv.model.out.plot)

SV area model with genotypes

dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15945 -0.04653  0.01473  0.03951  0.18609 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.015e-02  1.827e-02  -2.745   0.0092 ** 
## sv_area      4.426e-06  1.126e-07  39.315   <2e-16 ***
## group        1.813e-03  2.108e-02   0.086   0.9319    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06676 on 38 degrees of freedom
## Multiple R-squared:  0.9765, Adjusted R-squared:  0.9752 
## F-statistic: 788.8 on 2 and 38 DF,  p-value: < 2.2e-16
dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight, 
                                        group=genotype, color=genotype)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Genotype")
print(dry.sv.model.gen.plot)

Analyze VIS data

In this section we analyze geometric traits from the complete VIS data set.

Remove plants that were sampled for biomass measurements.

vis.data.zoom = vis.data.zoom[!vis.data.zoom$plant_id %in% st.data$plant_id,]

Predict fresh- and dry-weight biomass from linear models

vis.data.zoom$fw_biomass = predict.lm(object = sv.model, newdata=vis.data.zoom)
vis.data.zoom$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data.zoom)

Plant height analysis

Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.

height.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
                               vis.data.zoom$genotype == 'B100') &
                              (vis.data.zoom$treatment == 100 |
                               vis.data.zoom$treatment == 33),],
                     aes(x=dap, y=height_above_bound, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,60),
                                        name="Estimated height (cm)") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)

Plot height for all genotypes.

height.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
                                    vis.data.zoom$treatment == 33,],
                           aes(x=dap, y=height_above_bound, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,60),
                                        name="Estimated height (cm)") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)

Genotype differences in height

For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.

height_per_day = function(vis.data) {
  dap = c()
  genotype = c()
  int.low = c()
  int.up = c()
  est = c()
  
  # Loop through each day and genotype
  for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
    if(d != 24) {
      for(g in levels(factor(vis.data$genotype))) {
        h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
                          as.integer(vis.data$dap) == d,]$height_above_bound
        if(sd(h.data) != 0) {
          dap = c(dap,d)
          genotype = c(genotype,g)
          test = t.test(h.data)
          int.low = c(int.low,test$conf.int[1])
          int.up = c(int.up,test$conf.int[2])
          est = c(est,test$estimate)
        }
      }
    }
  }
  #drought_resp = drought_resp[-14,]
  results = data.frame(dap=as.numeric(dap),
                       genotype=genotype,
                       conf.int.low=int.low,
                       conf.int.up=int.up,
                       mean=est)  

  return(results)
}

Height per day 95% CI. Write results to file.

height.perday.results = height_per_day(vis.data.zoom)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)
##     dap genotype conf.int.low conf.int.up      mean
## 1    11      A10    4.1138046    5.278839  4.696322
## 2    11     B100    2.5259319    3.326661  2.926297
## 3    11     R102    1.2684481    3.282102  2.275275
## 4    11     R128    1.6661243    2.966349  2.316237
## 5    11     R133    2.8262168    5.442689  4.134453
## 6    11     R161    2.8964352    3.279606  3.088021
## 7    11     R187    2.2374378    4.006971  3.122204
## 8    11      R20    4.2648901    6.073519  5.169204
## 9    11      R70    3.3847556    4.309511  3.847133
## 10   11      R98    5.5427726    6.565899  6.054336
## 11   12      A10    4.6895754    6.770786  5.730181
## 12   12     B100    2.7955379    4.609936  3.702737
## 13   12     R102    2.0129905    5.055558  3.534274
## 14   12     R128    2.3173878    4.551756  3.434572
## 15   12     R133    0.1833279   16.040642  8.111985
## 16   12     R161    2.7075759    5.056310  3.881943
## 17   12     R187    4.3345350    6.973768  5.654151
## 18   12      R20  -12.0782382   28.023729  7.972745
## 19   12      R70    3.0529634    5.786161  4.419562
## 20   12      R98    7.0338485    9.288841  8.161345
## 21   13      A10    7.2633774    8.933993  8.098685
## 22   13     B100    5.4971969    6.483222  5.990209
## 23   13     R102    3.9400526    6.147069  5.043561
## 24   13     R128    4.8671575    6.363675  5.615416
## 25   13     R133    7.5927485    9.939040  8.765894
## 26   13     R161    5.6093371    7.035664  6.322501
## 27   13     R187    5.3732471    7.338795  6.356021
## 28   13      R20    7.2176543   10.572329  8.894992
## 29   13      R70    6.2027656    7.171094  6.686930
## 30   13      R98    8.7080117   11.816909 10.262460
## 31   14      A10    7.2992964   10.650554  8.974925
## 32   14     B100    7.2095969    8.594076  7.901836
## 33   14     R102    4.6619405    8.705043  6.683492
## 34   14     R128    5.6149445    9.147871  7.381408
## 35   14     R133   11.8870482   12.800615 12.343831
## 36   14     R161    5.5132820    8.133326  6.823304
## 37   14     R187    7.7707801   10.776425  9.273602
## 38   14      R20    4.4397218   17.446651 10.943186
## 39   14      R70    6.2590147   10.235183  8.247099
## 40   14      R98   11.0626447   12.903724 11.983184
## 41   15      A10   10.5250697   12.024864 11.274967
## 42   15     B100    8.6034387   10.159156  9.381297
## 43   15     R102    5.5131064    8.935807  7.224457
## 44   15     R128    8.3515328   10.617965  9.484749
## 45   15     R133   13.4816282   16.926665 15.204147
## 46   15     R161    8.5636321   10.277182  9.420407
## 47   15     R187    7.8971054   11.444819  9.670962
## 48   15      R20   12.5540041   16.579996 14.567000
## 49   15      R70    9.4534249   11.399830 10.426628
## 50   15      R98   14.4840908   15.932854 15.208472
## 51   16      A10   10.3183950   13.664952 11.991674
## 52   16     B100   10.3217123   12.069260 11.195486
## 53   16     R102    7.4126348   10.542533  8.977584
## 54   16     R128    8.2558336   15.124237 11.690035
## 55   16     R133   16.4716388   17.833934 17.152786
## 56   16     R161    8.2040259   11.746160  9.975093
## 57   16     R187   11.1664495   13.991223 12.578836
## 58   16      R20    4.4901376   25.384677 14.937407
## 59   16      R70   10.1166241   15.935022 13.025823
## 60   16      R98   16.3400141   17.763702 17.051858
## 61   17      A10   13.8288699   16.273493 15.051181
## 62   17     B100   12.8021984   14.751262 13.776730
## 63   17     R102    4.6784026   20.775747 12.727075
## 64   17     R128   13.3720766   17.669471 15.520774
## 65   17     R133   17.4574601   17.473869 17.465665
## 66   17     R161    9.2299092   20.579302 14.904606
## 67   17     R187   13.3249948   15.732300 14.528647
## 68   17      R20   17.2492825   17.540747 17.395015
## 69   17      R70   13.8775205   16.001331 14.939426
## 70   17      R98   17.4389051   17.482331 17.460618
## 71   18      A10   14.4288846   16.848835 15.638860
## 72   18     B100   15.1035778   16.693924 15.898751
## 73   18     R102   11.9392515   17.839681 14.889466
## 74   18     R128   14.4065180   16.592639 15.499579
## 75   18     R161   11.3654767   16.590256 13.977866
## 76   18     R187   14.2918010   18.966135 16.628968
## 77   18      R20   17.3669600   17.559323 17.463141
## 78   18      R70   16.2359873   18.115003 17.175495
## 79   18      R98   17.4695785   17.486983 17.478281
## 80   19      A10   15.5274345   17.339111 16.433273
## 81   19     B100   14.1045628   16.056801 15.080682
## 82   19     R102    9.2740593   18.382349 13.828204
## 83   19     R128   17.4682044   17.481868 17.475036
## 84   19     R133   17.4691725   17.478305 17.473739
## 85   19     R161   13.7376679   17.711630 15.724649
## 86   19     R187   15.1000979   17.975410 16.537754
## 87   19      R20   17.4722809   17.482904 17.477592
## 88   19      R70   16.3674982   17.915366 17.141432
## 89   20      A10   15.4815153   17.904691 16.693103
## 90   20     B100   14.9940006   18.062583 16.528292
## 91   20     R102   14.8917746   18.071450 16.481612
## 92   20     R128   15.8643804   18.103082 16.983731
## 93   20     R133   17.4653321   17.482145 17.473739
## 94   20     R161   13.7669020   17.208033 15.487467
## 95   20     R187   16.8032817   17.759659 17.281470
## 96   20      R20   17.4540444   17.497470 17.475757
## 97   20      R70   15.5391179   18.384945 16.962032
## 98   20      R98   17.4675530   17.483962 17.475757
## 99   21      A10   18.0979070   22.866964 20.482435
## 100  21     B100   17.2659492   20.740679 19.003314
## 101  21     R102   10.3649359   27.239881 18.802408
## 102  21     R128   23.2259512   29.985433 26.605692
## 103  21     R133   35.2919365   38.324288 36.808112
## 104  21     R161   15.6427963   23.752576 19.697686
## 105  21     R187   20.9302871   28.030502 24.480395
## 106  21      R20   28.5799826   34.970515 31.775249
## 107  21      R70   22.9352632   26.500120 24.717692
## 108  21      R98   28.6623854   35.443911 32.053148
## 109  22      A10   20.8063040   25.919929 23.363117
## 110  22     B100   17.8425711   25.156676 21.499623
## 111  22     R102   16.7256650   28.193792 22.459728
## 112  22     R128   22.0557320   29.847211 25.951471
## 113  22     R133   33.3894751   49.007692 41.198583
## 114  22     R161   16.4579791   20.578905 18.518442
## 115  22     R187   25.3813857   35.138840 30.260113
## 116  22      R20   23.4710805   58.736805 41.103943
## 117  22      R70   19.0935859   36.220406 27.656996
## 118  22      R98   33.3321134   40.188754 36.760433
## 119  23      A10   23.7887258   29.592135 26.690430
## 120  23     B100   20.1451618   25.383980 22.764571
## 121  23     R102   11.4174069   32.037610 21.727508
## 122  23     R128   29.2037138   35.281459 32.242586
## 123  23     R133   41.3614576   45.130164 43.245811
## 124  23     R161   20.4242532   25.227423 22.825838
## 125  23     R187   23.1617592   38.261040 30.711400
## 126  23      R20   34.0321909   44.150975 39.091583
## 127  23      R70   30.2497912   34.894308 32.572050
## 128  23      R98   31.3077114   41.034772 36.171242
## 129  25      A10   28.8418538   34.316538 31.579196
## 130  25     B100   21.7374532   31.214009 26.475731
## 131  25     R102   19.3852591   38.594612 28.989935
## 132  25     R128   30.4325931   40.029886 35.231239
## 133  25     R133   47.2030660   54.937143 51.070105
## 134  25     R161   21.5395173   26.577798 24.058658
## 135  25     R187   35.8052978   41.428971 38.617134
## 136  25      R20   24.4029064   70.696036 47.549471
## 137  25      R70   28.2682884   48.061912 38.165100
## 138  25      R98   40.7567043   49.171868 44.964286
## 139  26      A10   31.3246339   38.418576 34.871605
## 140  26     B100   26.7993073   34.574345 30.686826
## 141  26     R102   14.9726330   46.784895 30.878764
## 142  26     R128   35.9731220   44.331988 40.152555
## 143  26     R133   52.6693551   52.719093 52.694224
## 144  26     R161   23.0900303   39.985494 31.537762
## 145  26     R187   33.1394226   43.142960 38.141191
## 146  26      R20   49.2489502   53.411285 51.330117
## 147  26      R70   41.7251538   46.738979 44.232066
## 148  26      R98   34.7437085   49.954721 42.349215
## 149  27      A10   35.1602826   41.245993 38.203138
## 150  27     B100   27.0240390   37.926865 32.475452
## 151  27     R102   23.9486786   47.091606 35.520142
## 152  27     R128   35.9592946   47.623371 41.791333
## 153  27     R133   52.6385536   52.932629 52.785592
## 154  27     R161   24.2442800   30.727013 27.485647
## 155  27     R187   41.2496516   44.943105 43.096378
## 156  27      R20   48.8256683   55.318931 52.072300
## 157  27      R70   34.9404743   52.806382 43.873428
## 158  27      R98   47.5987571   53.069042 50.333900
## 159  28      A10   37.3444960   43.747626 40.546061
## 160  28     B100   31.0852941   39.304128 35.194711
## 161  28     R102   19.7785402   51.925226 35.851883
## 162  28     R128   45.7837167   50.042465 47.913091
## 163  28     R133   52.6330116   52.823749 52.728380
## 164  28     R161   28.0109976   44.448898 36.229948
## 165  28     R187   35.2088327   47.003038 41.105935
## 166  28      R20   52.5767078   52.970566 52.773637
## 167  28      R70   49.9879362   52.206923 51.097429
## 168  28      R98   42.0724187   54.742865 48.407642
## 169  29      A10   40.2283979   45.252515 42.740457
## 170  29     B100   33.4612448   46.507641 39.984443
## 171  29     R102   26.5575297   53.548336 40.052933
## 172  29     R128   40.3678421   53.396169 46.882005
## 173  29     R133   52.6874775   53.051070 52.869274
## 174  29     R161   31.3730339   37.445685 34.409360
## 175  29     R187   44.0292757   47.520641 45.774958
## 176  29      R20   52.6909028   52.987872 52.839387
## 177  29      R70   42.0424398   56.607023 49.324731
## 178  29      R98   52.0544743   52.857213 52.455843
## 179  30      A10   41.1590692   48.054219 44.606644
## 180  30     B100   36.0961707   44.787331 40.441751
## 181  30     R102   27.1510111   54.451174 40.801093
## 182  30     R128   51.7339576   52.709507 52.221732
## 183  30     R133   52.6812412   52.877987 52.779614
## 184  30     R161   37.2649692   46.088567 41.676768
## 185  30     R187   39.0026350   49.389770 44.196202
## 186  30      R20   52.7731212   52.959449 52.866285
## 187  30      R70   52.6255605   52.942207 52.783884
## 188  30      R98   48.3607663   54.312277 51.336522
## 189  31      A10   44.9409479   48.995137 46.968042
## 190  31     B100   39.6784848   50.930006 45.304246
## 191  31     R102   32.3953484   54.156046 43.275697
## 192  31     R128   45.9896516   54.628337 50.308994
## 193  31     R133   52.8090833   53.001192 52.905138
## 194  31     R161   35.3843979   41.461346 38.422872
## 195  31     R187   46.3300941   49.725217 48.027656
## 196  31      R20   52.5361938   53.321900 52.929047
## 197  31      R70   45.1991354   56.510708 50.854922
## 198  31      R98   52.7570164   52.941683 52.849349
## 199  32      A10   43.9442034   51.699444 47.821823
## 200  32     B100   42.0255866   49.295184 45.660385
## 201  32     R102   33.9379036   55.911966 44.924935
## 202  32     R128   53.2778901   54.523848 53.900869
## 203  32     R133   54.3529911   54.469816 54.411403
## 204  32     R161   39.2900732   49.085036 44.187554
## 205  32     R187   43.5507291   49.570772 46.560751
## 206  32      R20   54.4236761   54.468606 54.446141
## 207  32      R70   54.3062174   54.417340 54.361779
## 208  32      R98   53.3227703   54.787243 54.055007
## 209  33      A10   48.7441481   52.620165 50.682157
## 210  33     B100   45.2839248   53.441351 49.362638
## 211  33     R102   38.4873097   55.935265 47.211287
## 212  33     R128   48.7761391   55.972918 52.374528
## 213  33     R133   54.3633798   54.459427 54.411403
## 214  33     R161   41.5367102   45.138371 43.337541
## 215  33     R187   47.6165089   50.414019 49.015264
## 216  33      R20   54.3071074   54.578858 54.442983
## 217  33      R70   46.8503717   57.968160 52.409266
## 218  33      R98   54.3418212   54.470459 54.406140

Treatment differences in height

Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.

analyze_height = function(vis.data, genotype) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(vis.data$dap)))) {
    day = as.integer(day)
    control = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 100,]$height_above_bound
    drought = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 33,]$height_above_bound
    test = t.test(x=control, y=drought)
    days = c(days, day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Test for treatment effect on two-day intervals

a10.height.results = analyze_height(vis.data.zoom, 'A10')
b100.height.results = analyze_height(vis.data.zoom, 'B100')

Control for multiple testing by controlling the FDR

qvalues.height = p.adjust(c(a10.height.results$pvalue, b100.height.results$pvalue),method="fdr")
a10.height.results$qvalue = qvalues.height[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height[
  (nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]

Assign genotypes for merged table

a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11   0.77273960    2.643169 9.582125e-04 0.0032431807      A10
## 2   12  -0.29914242    2.422404 1.206740e-01 0.2097086778      A10
## 3   13  -0.08842057    2.115835 7.039538e-02 0.1534676662      A10
## 4   14  -0.32917060    2.580711 1.239188e-01 0.2097086778      A10
## 5   15   0.28426142    2.919397 1.924394e-02 0.0529208463      A10
## 6   16   0.14666528    4.378048 3.748556e-02 0.0868086644      A10
## 7   17  -0.87910116    3.018279 2.505424e-01 0.3674622138      A10
## 8   18   0.48991495    3.897893 1.460564e-02 0.0428432116      A10
## 9   19   1.15428894    4.393747 1.873909e-03 0.0058894273      A10
## 10  20   2.29234606    6.796709 2.596630e-04 0.0009520977      A10
## 11  21   3.42568864    8.778005 6.793406e-05 0.0003736373      A10
## 12  22   5.05021107   10.781276 3.072841e-06 0.0001352050      A10
## 13  23   4.77919478   12.771985 2.155551e-04 0.0008622202      A10
## 14  25   3.75450813   10.565738 1.699507e-04 0.0008308701      A10
## 15  26   4.41239580   12.298773 1.911306e-04 0.0008409746      A10
## 16  27   5.94525469   13.501681 1.946187e-05 0.0001951542      A10
## 17  28   5.54257526   13.398589 5.606520e-05 0.0003736373      A10
## 18  29   5.90446577   12.852781 7.632409e-06 0.0001679130      A10
## 19  30   6.17854616   13.985667 2.217662e-05 0.0001951542      A10
## 20  31   6.60003839   14.577970 1.300773e-05 0.0001907800      A10
## 21  32   6.21299029   15.161443 6.320013e-05 0.0003736373      A10
## 22  33   1.73822539   19.572166 2.690657e-02 0.0696405324      A10
print(b100.height.results)
##    dap conf.int.low conf.int.up     pvalue     qvalue genotype
## 1   11   -1.2594720  0.06992054 0.07673383 0.15346767     B100
## 2   12   -0.9886257  1.13640045 0.88658293 0.88658293     B100
## 3   13   -1.8408812  0.31467086 0.15621418 0.24547943     B100
## 4   14   -1.1123323  1.64762375 0.68794161 0.70394026     B100
## 5   15   -2.6455439  0.31836058 0.11604493 0.20970868     B100
## 6   16   -0.7864901  2.70009313 0.26085877 0.37025116     B100
## 7   17   -2.5824550  1.44347769 0.54836630 0.60320293     B100
## 8   18   -1.0746677  2.55761958 0.40054959 0.51835830     B100
## 9   19   -1.3873902  2.36907573 0.58999371 0.63316398     B100
## 10  20   -0.1953884  3.88187425 0.07458440 0.15346767     B100
## 11  21   -4.7055442  2.44925827 0.51436328 0.60320293     B100
## 12  22   -4.1699221  2.69196199 0.66036184 0.69180764     B100
## 13  23   -5.6805044  2.80840878 0.46303807 0.56593541     B100
## 14  25   -6.9139710  1.33553502 0.17774764 0.26968607     B100
## 15  26   -8.9304594 -0.52709508 0.02867192 0.07008691     B100
## 16  27   -7.3841917  1.12487918 0.14381616 0.23436707     B100
## 17  28   -8.7690652  0.63533080 0.08754828 0.16748367     B100
## 18  29   -6.1716464  2.70162092 0.43167803 0.54268095     B100
## 19  30   -6.3500065  2.33465144 0.35381611 0.48649716     B100
## 20  31   -2.6645326  4.95701544 0.54451691 0.60320293     B100
## 21  32   -2.8011143  5.16714436 0.54834969 0.60320293     B100
## 22  33   -3.1147374  7.53586129 0.38830762 0.51774350     B100

Output the A10 and B100 tables to the same file

write.table(a10.height.results, file='height.stats.csv', sep = ',',
            row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
            row.names = FALSE, append = TRUE, col.names = FALSE)

Biomass analysis

Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.

biomass.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
                                vis.data.zoom$genotype == 'B100') &
                               (vis.data.zoom$treatment == 100 |
                                vis.data.zoom$treatment == 33),],
                     aes(x=dap, y=fw_biomass, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(-1,41),
                                        name="Estimated biomass (g)") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)

Plot fresh-weight biomass for all genotypes.

biomass.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
                                    vis.data.zoom$treatment == 33,],
                           aes(x=dap, y=fw_biomass, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(-1,41),
                                        name="Estimated biomass (g)") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)

Treatment differences in biomass

Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.

analyze_biomass = function(vis.data, genotype) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(vis.data$dap)))) {
    day = as.integer(day)
    control = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 100,]$fw_biomass
    drought = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 33,]$fw_biomass
    # If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
    if (genotype == 'A10' | genotype == 'B100') {
      control = sample(control, size = length(drought))
    }
    test = t.test(x=control, y=drought)
    days = c(days, day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Test for treatment effect on two-day intervals

a10.biomass.results = analyze_biomass(vis.data.zoom, 'A10')
b100.biomass.results = analyze_biomass(vis.data.zoom, 'B100')

Control for multiple testing by controlling the FDR

qvalues.biomass = p.adjust(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue),method="fdr")
a10.biomass.results$qvalue = qvalues.biomass[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass[
  (nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) + 
                                 nrow(b100.biomass.results))]

Assign genotypes for merged table

a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11   0.06618061   0.2366926 1.407098e-03 2.948206e-03      A10
## 2   12  -0.08602355   0.2214458 3.733164e-01 4.211774e-01      A10
## 3   13  -0.14878889   0.2530395 6.000685e-01 6.286432e-01      A10
## 4   14  -0.11355314   0.4706694 2.189056e-01 2.832896e-01      A10
## 5   15  -0.16532964   0.7004685 2.138495e-01 2.832896e-01      A10
## 6   16  -0.13466603   1.1361927 1.146885e-01 1.802249e-01      A10
## 7   17  -0.71306484   2.5798564 2.404504e-01 3.022805e-01      A10
## 8   18   1.01265567   2.7833441 1.984776e-04 5.822009e-04      A10
## 9   19   0.59986296   2.8235229 3.885584e-03 7.433292e-03      A10
## 10  20   1.51836098   4.7372845 6.640184e-04 1.623156e-03      A10
## 11  21   2.11154234   6.0531431 3.243150e-04 8.918661e-04      A10
## 12  22   4.28169719   8.8377380 9.412210e-06 3.185671e-05      A10
## 13  23   4.94273400  10.4089661 6.069847e-05 1.907666e-04      A10
## 14  25   6.47172203  11.7791467 6.402658e-07 3.130189e-06      A10
## 15  26   7.10409340  11.8459095 3.113334e-08 1.712334e-07      A10
## 16  27   8.85914954  12.3725610 6.213467e-12 4.556542e-11      A10
## 17  28   9.61069065  12.9472564 5.285230e-13 5.813753e-12      A10
## 18  29   8.91289326  12.3369483 3.464972e-12 3.049175e-11      A10
## 19  30  10.04685735  12.9896112 3.742783e-14 1.646825e-12      A10
## 20  31  10.14299516  13.5602211 3.490344e-13 5.119172e-12      A10
## 21  32  10.04460532  13.3093630 1.626586e-13 3.578490e-12      A10
## 22  33   9.97305515  14.9878564 1.469255e-06 6.464720e-06      A10
print(b100.biomass.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11  -0.10273614  0.01603354 1.433897e-01 2.175569e-01     B100
## 2   12  -0.05331944  0.14763773 3.417760e-01 3.957407e-01     B100
## 3   13  -0.19433167  0.11252372 5.892521e-01 6.286432e-01     B100
## 4   14  -0.17223071  0.31125464 5.550045e-01 6.105049e-01     B100
## 5   15  -0.56141612  0.03872965 8.427560e-02 1.426203e-01     B100
## 6   16  -0.16839758  1.01001512 1.508052e-01 2.199402e-01     B100
## 7   17  -1.20761685  0.83516601 7.032150e-01 7.195688e-01     B100
## 8   18  -0.16852448  1.69244380 1.035239e-01 1.687056e-01     B100
## 9   19  -1.11131797  1.07162908 9.703722e-01 9.703722e-01     B100
## 10  20   0.95683425  3.18300525 8.908319e-04 2.062979e-03     B100
## 11  21  -1.02609382  3.09523737 3.076201e-01 3.759801e-01     B100
## 12  22  -0.71744932  4.15312089 1.549578e-01 2.199402e-01     B100
## 13  23  -1.72716481  7.54585618 1.698545e-01 2.335500e-01     B100
## 14  25   0.30686251  6.13591628 3.222730e-02 5.672006e-02     B100
## 15  26  -2.16265890  6.18863717 3.194544e-01 3.798917e-01     B100
## 16  27   4.43237220  9.24389988 8.819553e-06 3.185671e-05     B100
## 17  28   1.33865042  8.29025819 9.646199e-03 1.768470e-02     B100
## 18  29   2.22560476  8.60958217 2.396530e-03 4.793060e-03     B100
## 19  30   2.93818406 10.03007759 1.320144e-03 2.904316e-03     B100
## 20  31   8.43451895 12.61765594 5.634832e-10 3.541894e-09     B100
## 21  32   4.53402087 12.71198560 4.007013e-04 1.037109e-03     B100
## 22  33   8.36684530 14.24962413 1.865143e-06 7.460572e-06     B100

Output the A10 and B100 tables to the same file

write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
            row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
            row.names = FALSE, append = TRUE, col.names = FALSE)

Biomass response to drought

Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.

drought_response = function(vis.data, genotypes) {
  # Initialize drought response data frame with days
  drought_resp = data.frame(day=c(
    min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
  
  # Initialize genotypes
  for(g in genotypes) {
    control = paste(g,'control',sep='')
    drought = paste(g,'drought',sep='')
    
    # Calculate median biomass per treatment per day
    drought_resp[,paste(control)] = 0
    drought_resp[,paste(drought)] = 0
    for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
      drought_resp[drought_resp$day == d, paste(control)] = 
        median(vis.data[vis.data$genotype == g &
                        vis.data$treatment == 100 &
                        as.integer(vis.data$dap) == d,]$fw_biomass)
      drought_resp[drought_resp$day == d, paste(drought)] = 
        median(vis.data[vis.data$genotype == g &
                        vis.data$treatment == 33 &
                        as.integer(vis.data$dap) == d,]$fw_biomass)
    }
  }
  drought_resp = drought_resp[-14,]
  
  # Calculate biomass loss to drought
  days = c()
  response = c()
  genotype=c()
  for(r in 1:nrow(drought_resp)) {
    days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
    for(g in genotypes) {
      control = paste(g,'control',sep='')
      drought = paste(g,'drought',sep='')
      response = c(response, drought_resp[r,paste(drought)] -
                   drought_resp[r,paste(control)])
      genotype = c(genotype, g)
    }
  }
  dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
  return(dr.df)
}

A10 vs B100 drought responses

drought.response.parents = drought_response(vis.data.zoom,c('A10','B100'))

Plot drought response results

genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
                   aes(x=day, y=response, group=genotype, color=genotype)) +
                   geom_point(size=2.5) +
                   geom_smooth(method="loess") +
                   scale_x_continuous(name="Days after planting") +
                   scale_y_continuous(name="Reduced accumulated biomass (g)") +
                   theme_bw() +
                   theme(legend.position=c(0.8,0.8),
                         axis.title.x=element_text(face="bold"),
                         axis.title.y=element_text(face="bold")) +
                   scale_colour_manual(values=genotype.colors) +
                   labs(color='Genotype')
resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)

Growth curve analysis

In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.

First load three functions for non-linear growth curve analysis reproduced from:

Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.

output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
  coef <- coef(fit)
  params <- transform_param.logis(coef)
  rates <- list()
  groups <- rownames(params)
  n.groups <- nrow(coef)
  
  # compute rates for each group seperately
  rates <- list()
  for(i in 1:(n.groups)){
    K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
    rates[[i]] = data.frame(
      times = times,
      M    = (M0*K)/(M0+(K-M0)*exp(-r*times)),
      AGR  = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
    )
    rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
    rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
    if(LOG == T){
      rates[[i]]$RGRt <- rates[[i]]$AGR
      rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
      rates[[i]]$AGR  <- rates[[i]]$AGR*exp(rates[[i]]$M)
    }
    # commute CIs for each group's estaimates, if desired
    if(CI == T){
      cov   <- summary(fit)$cov[i,,]
      x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
      x$K  <- y[,1]
      x$r  <- 1/y[,3]
      x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3])) 
      M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
      for(j in 1:nrow(x)){
        K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
        M[j,]     <- (M0*K)/(M0+(K-M0)*exp(-r*times))
        AGR[j,]   <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
        RGRt[j,]  <- AGR[j,]/M[j,]
        RGRm[j,]  <- r*(1 - M[j,]/K)
        if(LOG ==T){
          RGRt[j,] <- AGR[j,]
          RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
          AGR[j,]  <- AGR[j,]*exp(M[j,])
        }
      }
    }
    CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
    rates[[i]]  <-  cbind(rates[[i]], CIs)
  }
  names(rates) <- rownames(params)
  
  # now compute differences among groups
  diffs <- list()
  for(i in 1:(n.groups-1)){
    Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
    for(j in (i+1):n.groups){
      Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
      diffs.ij = data.frame(
        times = times,
        diffM     = rates[[i]]$M    - rates[[j]]$M,
        diffAGR   = rates[[i]]$AGR  - rates[[j]]$AGR,
        diffRGRt  = rates[[i]]$RGRt - rates[[j]]$RGRt
      )
      # comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
      Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
      Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
      #diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
      diffs.ij$Mseq <- seq(Mmin, Mmax, length = length(times))
      if(LOG == F){
        diffs.ij$diffRGRm  <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
      } else{
        diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
      }    
      
    }
    if(CI == T){
      # get params for group i
      covi   <- summary(fit)$cov[i,,]
      xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
      xi$K  <- yi[,1]
      xi$r  <- 1/yi[,3]
      xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3])) 
      
      # get params for group j
      covj   <- summary(fit)$cov[j,,]
      xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
      xj$K  <- yj[,1]
      xj$r  <- 1/yj[,3]
      xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3])) 
      
      # now compute diffs for each random set of drawn parameters
      Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
      for(k in 1:nrow(xi)){
        Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
        Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
        Mi[k,]    <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
        Mj[k,]    <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
        AGRi[k,]  <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
        AGRj[k,]  <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
        RGRti[k,] <- AGRi[k,]/Mi[k,]
        RGRtj[k,] <- AGRj[k,]/Mj[k,]
        RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
        RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
        if(LOG == T){
          RGRti[k,] <- AGRi[k,]
          RGRtj[k,] <- AGRj[k,]
          RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
          RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
          AGRi[k,]  <- AGRi[k,]*exp(Mi[k,])
          AGRj[k,]  <- AGRj[k,]*exp(Mj[k,])
        }
        diffM[k,]    <- Mi[k,]    - Mj[k,]
        diffAGR[k,]  <- AGRi[k,]  - AGRj[k,]
        diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
        diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
      }
      CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
      diffs[[paste(groups[i], groups[j], sep = "_")]]  <-  cbind(diffs.ij, CIs)
    } else{
      diffs[[paste(groups[i], groups[j], sep = "_")]]  <-  diffs.ij
    }
  } # end loop over pairwise combinations of groups
  
  out <- list(params = params, rates = rates, diffs = diffs)
  return(out)
}

transform_param.logis <- function(coef){
  K = coef[1]
  r = 1/(coef[3])
  M0 =  K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
  if(is.data.frame(K)){
    out <- cbind(K, r, M0)
  } else {
    out <- c(K, r, M0)
  }
  names(out) <- c("K", "r", "M0")
  return(out)
}

# this function returns confidence envelopes around growth trajectories, and growth rates. 
summarizer <- function(dat, alpha){
  n <- length(dat)
  quantiles <- c(alpha/2, 1-(alpha/2))
  CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
  names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
  for(i in 1:n){
    CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]],    2, quantile, quantiles, na.rm = T))
  }
  return(CIs)
}

Subset the VIS data

sub.vis.data = vis.data.zoom[(vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),
                        c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)

Group data

grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)

Fit three-component logistic functions for each genotype x treatment group

fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
                    data = grouped.sub.vis.data)

Output estimated growth rate parameters at even time intervals

est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
                                     CI = TRUE, LOG = FALSE, alpha = 0.05)

Compute the time and magnitude of maximum growth rate

logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
logis.results$AGR.lo = 0
logis.results$AGR.hi = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
  logis.results[logis.results == group.order[i],]$max.AGR.dap = 
    out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR == 
                                             max(out.fit.logis$rates[[i]]$AGR)]
  logis.results[logis.results == group.order[i],]$max.AGR =
    max(out.fit.logis$rates[[i]]$AGR)
  
  logis.results[logis.results == group.order[i],]$AGR.lo = 
    out.fit.logis$rates[[i]]$AGR.lo[out.fit.logis$rates[[i]]$AGR == 
                                             max(out.fit.logis$rates[[i]]$AGR)]
  
  logis.results[logis.results == group.order[i],]$AGR.hi = 
    out.fit.logis$rates[[i]]$AGR.hi[out.fit.logis$rates[[i]]$AGR == 
                                             max(out.fit.logis$rates[[i]]$AGR)]
}
print(logis.results)
##       group  max.AGR max.AGR.dap   AGR.lo   AGR.hi
## 1   A10-100 3.031821    24.23556 2.962031 3.099600
## 2    A10-33 1.800921    25.36054 1.724413 1.889975
## 3  B100-100 2.428901    26.48553 2.378344 2.479694
## 4   B100-33 1.693055    24.68555 1.604441 1.784634
## 5  R102-100 2.679992    28.06050 2.596513 2.757770
## 6   R102-33 1.890607    26.48553 1.813447 1.977527
## 7  R128-100 3.090714    25.58554 3.003102 3.174770
## 8   R128-33 1.847410    24.46055 1.770401 1.930073
## 9  R133-100 3.380220    23.11057 3.222844 3.526918
## 10  R133-33 2.292966    23.33557 2.180286 2.421985
## 11 R161-100 2.952080    25.58554 2.853000 3.043931
## 12  R161-33 1.813290    24.68555 1.729947 1.913824
## 13 R187-100 2.535501    24.46055 2.451728 2.629045
## 14  R187-33 1.676044    24.46055 1.596170 1.763532
## 15  R20-100 3.360632    22.88557 3.225079 3.499300
## 16   R20-33 1.949526    22.66058 1.858486 2.053219
## 17  R70-100 2.978971    24.01056 2.886229 3.083468
## 18   R70-33 1.654110    23.56057 1.571712 1.742591
## 19  R98-100 2.597831    24.23556 2.507745 2.691376
## 20   R98-33 1.538670    24.68555 1.467270 1.612429

Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time

parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")

# Biomass over time
gca.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
                                 vis.data.zoom$genotype == 'B100') &
                 (vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),],
                  aes(x=dap, y=fw_biomass, color=factor(group))) +
                  geom_point(size=2.5) +
                  scale_x_continuous(name="Days after planting") +
                  scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  scale_colour_manual(name='Genotype-treatment',
                                      values=group.colors, labels=parent.groups)

agr.plot = ggplot() + 
           scale_x_continuous(name="Days after planting") +
           scale_y_continuous(name="Absolute growth rate (g/day)") +
           theme_bw() +
           theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
           scale_colour_manual(name='Genotype-treatment',
                                      values=group.colors, labels=parent.groups)

# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
  group = parent.groups[g]
  i = match(group, group.order)
  group.rates = as.data.frame(out.fit.logis$rates[i])
  col.name = gsub('-', '.', group)
  conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
  agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
  gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
                                     fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.M', sep='')),
              color=group.colors[g])
  
  agr.plot = agr.plot +
    geom_polygon(data=agr.conf.int, aes(x=x, y=y),
                 fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.AGR', sep='')),
              color=group.colors[g])
}
gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)

agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)

Plot logistic growth models for all Setaria genotypes

plot_gca = function(genotype, vis.data, out.fit.logis) {
  groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
  group.colors = c("#00BFC4", "#F8766D")
  # Biomass over time
  gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
                   (vis.data$treatment == 100 | vis.data$treatment == 33),],
                    aes(x=dap, y=fw_biomass, color=factor(treatment))) +
                    geom_point(size=1) +
                    scale_x_continuous(name="Days after planting") +
                    scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
                    theme_bw() +
                    theme(legend.position='none',
                          axis.title.x=element_text(face="bold"),
                          axis.title.y=element_text(face="bold")) +
                    labs(color="Treatment") +
                    facet_grid(. ~ genotype)
  # Add logistic growth models and confidence intervals
  for(g in 1:length(groups)) {
    group = groups[g]
    i = match(group, group.order)
    group.rates = as.data.frame(out.fit.logis$rates[i])
    col.name = gsub('-', '.', group)
    conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
    gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
                                     fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.M', sep='')),
              color=group.colors[g])
  }
  return(gca.plot)
}

Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/

multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

Generate the Setaria biomass/growth plots

a10.gca.plot = plot_gca('A10', vis.data.zoom, out.fit.logis)
b100.gca.plot = plot_gca('B100', vis.data.zoom, out.fit.logis)
r102.gca.plot = plot_gca('R102', vis.data.zoom, out.fit.logis)
r128.gca.plot = plot_gca('R128', vis.data.zoom, out.fit.logis)
r133.gca.plot = plot_gca('R133', vis.data.zoom, out.fit.logis)
r161.gca.plot = plot_gca('R161', vis.data.zoom, out.fit.logis)
r187.gca.plot = plot_gca('R187', vis.data.zoom, out.fit.logis)
r20.gca.plot = plot_gca('R20', vis.data.zoom, out.fit.logis)
r70.gca.plot = plot_gca('R70', vis.data.zoom, out.fit.logis)
r98.gca.plot = plot_gca('R98', vis.data.zoom, out.fit.logis)

Supplemental Figure S3

multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot, 
          r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot, 
          r70.gca.plot, r98.gca.plot, cols=3)

agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)

Water use efficiency

In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.

Download water volume data (n = 33,496 water applications).

if (!file.exists('B2_watering_data_phenofront.csv')) {
  download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
                'B2_watering_data_phenofront.csv')
}

Read data and format for analysis

water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)

Remove the empty pot controls (Barcodes start with zero).

water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]

Remove snapshots that were not valid waterings (water.amount = -1).

water.data = water.data[water.data$water.amount != -1,]

Add a treatment group column coded in the barcode.

water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66

Add a plant genotype column coded in the barcode.

water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'

Add a genotype x treatment group column.

water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')

Encode the calendar time column as a date-time.

water.data$timestamp = ymd_hms(water.data$timestamp)

Add a column for days after planting since the planting date.

water.data$dap = as.numeric(water.data$timestamp - planting_date)

Some plants were removed from the system but were not inactivated, so remove them.

bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
                                     water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]

Data for Figure 1B

Extract water data for S. viridis for Figure 1B.

sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
                       water.data$genotype == 'A10',]

Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.

sv.water$early.late <- NA  #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
                    hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]

Plot the S. viridis water volume data

water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
                                  group = factor(interaction(early.late,
                                                             treatment,dap)),
                                  fill = factor(interaction(early.late,
                                                            treatment)),
                                  color = factor(interaction(early.late,
                                                             treatment)))) +
                     geom_boxplot(outlier.colour = NULL) +
                     scale_x_continuous("Days after planting") +
                     scale_y_continuous("Water volume (ml)") +
                     theme_bw() +
                     theme(legend.position = c(0.2,0.8),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(fill='Water treatment', color='Water treatment')
water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)

WUE function

Calculate WUE for S. viridis and S. italica

wue.parents = function(water, biomass) {
  # WUE data vectors
  dap.list = c()
  plant.list = c()
  water.list = c()
  biomass.list = c()
  treatment.list = c()
  group.list=c()
  
  # Get unique barcodes for biomass
  barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
                            (biomass$treatment == 100 |
                             biomass$treatment == 33),]$plant_id)
  for(barcode in barcodes) {
    snapshots = biomass[biomass$plant_id==barcode,]
    snapshots = snapshots[with(snapshots, order(dap)),]
    for(row in 1:nrow(snapshots)) {
      total.water = sum(water[water$dap <= snapshots[row,]$dap &
                              water$plant.barcode == barcode,]$water.amount)
      
      dap.list = c(dap.list, snapshots[row,]$dap)
      plant.list = c(plant.list, barcode)
      water.list = c(water.list, total.water)
      biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
      treatment.list = c(treatment.list, snapshots[row,]$treatment)
      group.list = c(group.list,snapshots[row,]$group)
    }
  }
  
  wue.data = data.frame(plant_id=plant.list,
                        dap=dap.list,
                        water=water.list,
                        biomass=biomass.list,
                        treatment=treatment.list,
                        group=group.list)
  return(wue.data)
}

parents.wue = wue.parents(water.data, vis.data.zoom)
parents.wue = parents.wue[parents.wue$water != 0,]

Plot S. viridis and S. italica WUE in mg/mL.

wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
                                   color=factor(group))) +
                  geom_point(size=2.5) +
                  geom_smooth(method="loess", size=1) +
                  scale_x_continuous(name="Days after planting") +
                  scale_y_continuous(lim=c(-2.1, 26),
                                     name="Water-use efficiency (mg/mL)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  labs(color='Genotype-treatment')
wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)

Statistical analysis of WUE

Genotype WUE function

wue = function(water, biomass, genotype) {
  # WUE data vectors
  dap.list = c()
  plant.list = c()
  water.list = c()
  biomass.list = c()
  treatment.list = c()
  
  # Get unique barcodes for biomass
  barcodes = unique(biomass[biomass$genotype==genotype &
                           (biomass$treatment == 100 |
                            biomass$treatment == 33),]$plant_id)
  for(barcode in barcodes) {
    snapshots = biomass[biomass$plant_id==barcode,]
    snapshots = snapshots[with(snapshots, order(dap)),]
    for(row in 1:nrow(snapshots)) {
      total.water = sum(water[water$dap <= snapshots[row,]$dap &
                              water$plant.barcode == barcode,]$water.amount)
      
      dap.list = c(dap.list, snapshots[row,]$dap)
      plant.list = c(plant.list, barcode)
      water.list = c(water.list, total.water)
      biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
      treatment.list = c(treatment.list, snapshots[row,]$treatment)
    }
  }
  
  wue.data = data.frame(plant_id=plant.list,
                        dap=dap.list,
                        water=water.list,
                        biomass=biomass.list,
                        treatment=treatment.list)
  return(wue.data)
}

Analyse pair-wise treatment differences for two-day increments for each genotype.

analyze_wue = function(wue.data) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(wue.data$dap)))) {
    day = as.integer(day)
    control = wue.data[(as.integer(wue.data$dap) == day |
                        as.integer(wue.data$dap) == day + 1) &
                        wue.data$treatment == 100,]
    drought = wue.data[(as.integer(wue.data$dap) == day |
                        as.integer(wue.data$dap) == day + 1) &
                        wue.data$treatment == 33,]
    control.wue = control$biomass / control$water
    drought.wue = drought$biomass / drought$water
    test = t.test(x=control.wue,y=drought.wue)
    days = c(days,day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Analyze WUE for each genotype and aggregate results

wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
  genotype.wue = wue(water.data, vis.data.zoom, genotype)
  genotype.wue = genotype.wue[genotype.wue$water != 0,]
  genotype.wue.results = analyze_wue(genotype.wue)
  genotype.wue.results$genotype = genotype
  wue.results = rbind(wue.results, genotype.wue.results)
}

Control for multiple testing by controlling the FDR

qvalues.wue = p.adjust(wue.results$pvalue,method="fdr")
wue.results$qvalue = qvalues.wue
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)
##    dap  conf.int.low   conf.int.up       pvalue genotype     qvalue
## 1   11  0.0001681857  7.003870e-04 0.0023589893      A10 0.03459851
## 2   12 -0.0002932779  6.855075e-04 0.4161121704      A10 0.79604067
## 3   13 -0.0003981046  7.310925e-04 0.5504252348      A10 0.85019455
## 4   14 -0.0003666556  1.086597e-03 0.3151858236      A10 0.65930472
## 5   15 -0.0006318583  1.290448e-03 0.4824013893      A10 0.85019455
## 6   16 -0.0011579704  1.505706e-03 0.7841660738      A10 0.85019455
## 7   17 -0.0030117768  2.033894e-03 0.6757645778      A10 0.85019455
## 8   18 -0.0019852804  1.537689e-03 0.7920861585      A10 0.85019455
## 9   19 -0.0022189870  1.423300e-03 0.6543284916      A10 0.85019455
## 10  20 -0.0018767476  2.487869e-03 0.7709120024      A10 0.85019455
## 11  21 -0.0011815956  3.333287e-03 0.3296523601      A10 0.65930472
## 12  22 -0.0011941099  3.447888e-03 0.3180544763      A10 0.65930472
## 13  23 -0.0030785137  4.736555e-03 0.6313454757      A10 0.85019455
## 14  25 -0.0017208810  2.529517e-03 0.6913194980      A10 0.85019455
## 15  26 -0.0016761719  2.370431e-03 0.7214693240      A10 0.85019455
## 16  27 -0.0015099451  1.887345e-03 0.8191576027      A10 0.85816511
## 17  28 -0.0016265463  1.712963e-03 0.9574064089      A10 0.95740641
## 18  29 -0.0016959640  1.418049e-03 0.8545599052      A10 0.87443339
## 19  30 -0.0017780029  1.240136e-03 0.7142953563      A10 0.85019455
## 20  31 -0.0018270439  1.050066e-03 0.5808623905      A10 0.85019455
## 21  32 -0.0018669052  9.419976e-04 0.5014284454      A10 0.85019455
## 22  33 -0.0019935110  2.547949e-03 0.7922267419      A10 0.85019455
## 23  11 -0.0004686284  1.653400e-04 0.3231119733     B100 0.65930472
## 24  12 -0.0003176706  4.659287e-04 0.6970305202     B100 0.85019455
## 25  13 -0.0006793007  2.052254e-04 0.2809421281     B100 0.65930472
## 26  14 -0.0005693316  7.964500e-04 0.7298612911     B100 0.85019455
## 27  15 -0.0013429980  2.693045e-04 0.1789099092     B100 0.56228829
## 28  16 -0.0011018089  1.726678e-03 0.6393283914     B100 0.85019455
## 29  17 -0.0053976734 -4.565038e-04 0.0249473645     B100 0.11089327
## 30  18 -0.0036398672  2.934236e-04 0.0894445265     B100 0.30273532
## 31  19 -0.0050384129 -6.422642e-04 0.0147704390     B100 0.08123741
## 32  20 -0.0044244661 -5.891914e-05 0.0447774175     B100 0.17910967
## 33  21 -0.0064021316 -1.321929e-03 0.0053896550     B100 0.03952414
## 34  22 -0.0057466915 -1.176396e-03 0.0051207769     B100 0.03952414
## 35  23 -0.0062955680 -1.663831e-03 0.0033048057     B100 0.03635286
## 36  25 -0.0058026119 -1.706030e-03 0.0009451266     B100 0.03381652
## 37  26 -0.0052016475 -1.378281e-03 0.0015371147     B100 0.03381652
## 38  27 -0.0043036325 -5.750897e-04 0.0125215176     B100 0.07870668
## 39  28 -0.0037893424 -2.728566e-04 0.0252030152     B100 0.11089327
## 40  29 -0.0033915507  1.520038e-04 0.0711932892     B100 0.26104206
## 41  30 -0.0027631015  6.920535e-04 0.2276453890     B100 0.65930472
## 42  31 -0.0025651583  7.463200e-04 0.2647434581     B100 0.65930472
## 43  32 -0.0025173692  8.428216e-04 0.3102535812     B100 0.65930472
## 44  33 -0.0029600383  2.256215e-03 0.7752593144     B100 0.85019455

Model tiller number

In this section we use linear modeling to estimate tiller number.

Calculate height-width ratio.

vis.data.zoom$height_width_ratio = vis.data.zoom$height_above_bound / vis.data.zoom$extent_x

Download data for manually measured tiller counts 195 random images.

if (!file.exists('200tiller_counts.txt')) {
  download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
                '200tiller_counts.txt')
}

Read manual tiller count data.

tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)

Format date.

tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
                                            tiller.counts$month,
                                            tiller.counts$day, sep='-'),
                                      paste(tiller.counts$hours,
                                            tiller.counts$minutes,
                                            tiller.counts$seconds, sep=':'), sep=' '))

Merge with VIS snapshot table.

tiller.counts = merge(x=tiller.counts, y=vis.data.zoom, by = 'date')

Model tiller counts.

tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
                         extent_x + height_width_ratio, data=tiller.counts)

Variance inflation factor.

vif(tiller.model.full)
##         fw_biomass height_above_bound           solidity 
##          16.615475          21.832350           1.591739 
##           extent_x height_width_ratio 
##          14.738418           3.238483

Drop height and recalculate VIF

tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
                     height_width_ratio, data=tiller.counts)
vif(tiller.model1)
##         fw_biomass           solidity           extent_x 
##          10.014163           1.363769          11.118321 
## height_width_ratio 
##           1.847957

Drop extent x.

tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
                   data=tiller.counts)
vif(tiller.model2)
##         fw_biomass           solidity height_width_ratio 
##           1.024829           1.041759           1.039990
summary(tiller.model2)
## 
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio, 
##     data = tiller.counts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6877 -1.2636 -0.0169  0.8454  5.6053 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.23276    0.54698   9.567  < 2e-16 ***
## fw_biomass          0.21998    0.01254  17.538  < 2e-16 ***
## solidity            0.15006    1.91657   0.078    0.938    
## height_width_ratio -2.19370    0.32640  -6.721 2.12e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.707 on 187 degrees of freedom
## Multiple R-squared:  0.644,  Adjusted R-squared:  0.6382 
## F-statistic: 112.7 on 3 and 187 DF,  p-value: < 2.2e-16

Drop solidity

tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)
## 
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6886 -1.2680 -0.0274  0.8433  5.6022 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.26031    0.41771  12.593  < 2e-16 ***
## fw_biomass          0.21986    0.01242  17.709  < 2e-16 ***
## height_width_ratio -2.18932    0.32072  -6.826 1.17e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.703 on 188 degrees of freedom
## Multiple R-squared:  0.6439, Adjusted R-squared:  0.6402 
## F-statistic:   170 on 2 and 188 DF,  p-value: < 2.2e-16

Plot partial regressions.

avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")

Download data for manually measured tiller counts for 646 random images.

if (!file.exists('660tiller_counts.txt')) {
  download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
                '660tiller_counts.txt')
}

Predict tiller count for a second set of ~660 randomly selected plants.

tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data.zoom, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
                                              newdata=tiller660)

Summarize the difference between the predicted and manual tiller counts

tiller.diff = tiller660$predicted_tiller_count - tiller660$tiller_count
t.test(tiller.diff)
## 
##  One Sample t-test
## 
## data:  tiller.diff
## t = 14.7453, df = 645, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.7978015 1.0429342
## sample estimates:
## mean of x 
## 0.9203678

Plot manual versus predicted tiller counts.

tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
                                         color=factor(genotype))) +
                     geom_point(size=2.5) +
                     geom_abline(intercept=0, slope=1) +
                     scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
                     scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.45,0.95),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold"),
                           legend.direction="horizontal") +
                     labs(color="Genotype")
tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)

Predict tiller counts for the whole data set.

vis.data.zoom$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data.zoom)

Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.

tc.par.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
                                vis.data.zoom$genotype == 'B100') &
                               (vis.data.zoom$treatment == 100 |
                                vis.data.zoom$treatment == 33),],
                     aes(x=dap, y=tiller_count, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,12),
                                        name="Estimated tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).

Plot estimated tiller counts for all genotypes.

tiller.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
                                    vis.data.zoom$treatment == 33,],
                           aes(x=dap, y=tiller_count, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,12),
                                        name="Estimated tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).

Trait table

Output a table of all of the traits input and added here. Export all of the VIS traits (except color).

h.table = vis.data.zoom
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
  row = h.table[r,]
  total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
  h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)